Simulating Strongly Correlated Quantum Systems with Tree Tensor Networks 
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We present a tree-tensor-network-based method to study strongly correlated systems with non- 
local interactions in higher dimensions. Although the momentum-space and quantum-chemistry 
versions of the density matrix renormalization group (DMRG) method have long been applied to 
such systems, the spatial topology of DMRG-based methods allows efficient optimizations to be 
carried out with respect to one spatial dimension only. Extending the matrix-product-state picture, 
we formulate a more general approach by allowing the local sites to be coupled to more than two 
neighboring auxiliary subspaces. Following Shi. et. al. [Phys. Rev. A, 74, 022320 (2006)], we treat 
a tree-like network ansatz with arbitrary coordination number z, where the z = 2 case corresponds 
to the one-dimensional scheme. For this ansatz, the long-range correlation deviates from the mean- 
field value polynomially with distance, in contrast to the matrix-product ansatz, which deviates 
exponentially. The computational cost of the tree-tensor-network method is significantly smaller 
than that of previous DMRG-based attempts, which renormalize several blocks into a single block. 
In addition, we investigate the effect of unitary transformations on the local basis states and present 
a method for optimizing such transformations. For the 1-d interacting spinless fermion model, the 
optimized transformation interpolates smoothly between real space and momentum space. Calcu- 
lations carried out on small quantum chemical systems support our approach. 



I. INTRODUCTION 

Understanding and simulating strongly correlated sys- 
tems has long been a major challenge in theoreti- 
cal physics and in theoretical chemistry. In the past 
two decades, the density-matrix renormalization group 
(DMRG) method has been applied effectively to study 
problems in these fields*^ 2 - In particular, it has been 
widely used to study fermionic and spin-chain prob- 
lems in one dimension for models with both local 
and long-range interactions. Application to systems 
with long-range interactions gained impetus when the 
method was reformulated to treat models defined in mo- 
mentum spaced - — (MS-DMRG) and quantum chemistry 
calculations^— (QC-DMRG). Common characteristics 
of these approaches are that a higher dimensional sys- 
tem is mapped to a one-dimensional chain and that vari- 
ational approximations to the cigenstates are obtained 
by an iterative diagonalization procedure. Introduction 
of various quantum information entropie s 5 ' 17 " — and the 
reformulation of the problem in terms of matrix product 
states22r— (MPS) has clarified the mathematical under- 
pinnings of the method. 

Recently, extensions to two dimensions based on the 
controlled manipulation of entanglement between sub- 
systems has led to alternate methods that can be 
viewed as generalizations of matrix-product-state-based 
methods^ These methods include projected entangled 
pair state o 20 ' 24 " — (PEPS), the multiscale-entanglement- 
renormalization ansata 2 ^ (MERA), and correlator prod- 
uct states^ (CPS) or complete-graph tensor network 
states^ (CGTN). The PEPS and MERA methods have 
already shown considerable promise for frustrated and 



fermionic problems; they do not suffer from the fermion 
sign problem that appears in quantum Monte Carlo 
simulations! 30 ' 31 These new methods, however, have pri- 
marily been restricted to the treatment of local Hamil- 
tonians. Thus, developing effective algorithms to treat 
higher dimensional systems in which the interactions are 
nonlocal remains an important problem. 

In order to efficiently treat quantum chemical systems, 
it has become evident in the past few years that methods 
must take into account more general spatial topology. It 
has been shown that quantum information entropy^ 5 - 




FIG. 1. (Color online) Tree tensor network with z = 3 (a) 
and 2 = 4 (b). The structure of the tensors is shown in (c) 
and (d). The bonds indicate the virtual indices ai,...,a z 
and the circle the physical index k. 
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can be used to determine the entanglement or quantum 
correlation among sites or orbitals in a pairwise way. 
Such a two-dimensional entanglement matrix leads to a 
picture of the topology of orbitals that corresponds to a 
multiply connected network. Two possible approaches to 
optimize computational efficiency are to work with an ap- 
propriate network-like structure that reflects the entan- 
glement topology and to vary the single-particle basis^ 
so that entanglement is reduced. For states based on 
a one-dimensional topology, the only possibility to take 
entanglement topology into account is to reorder the or- 
bitals. The next step in generalizing the topology is to 
use a tree network. This was first pointed out in Ref. [33[ . 
A tree network makes it possible to reduce the distance 
between highly entangled orbitals when they are multiply 
connected. In addition, the coordination number z can 
be varied from orbital to orbital to adapt the variational 
state to a particular entanglement topology. For a tree 
network that is bipartite, it is possible to adapt methods 
and optimizations developed for matrix-product-state- 
based algorithms such as the DMRG^i to tensor-product 
states on the tree network. This is because the tree net- 
work method can be related to a generalized DMRG with 
z blocks instead of two. In addition, since an effective 
Hamiltonian can be formed, real time evolution, like that 
done in t-DMRG-^^— and in PEPS^, can be carried 
out. Real time evolution could be used, for example, to 
calculate spectral functions for quantum chemical sys- 
tems. 

The use of more general spatial topologies is also po- 
tentially important for quantum impurity problems. In 
these systems, the impurity subsystem, which can con- 
sist of either a single site or a strongly correlated clus- 
ter, is coupled to a bath of free fermions. The classical 
method to treat quantum impurity problems numerically 
is the numerical renormalization group (NRG)2£. The 
starting point of the NRG method is the mapping of the 
problem onto a one-dimensional semi-infinite lattice in 
which the first site (or set of sites) describes the impu- 
rity subsystem, and the remainder of the chain represents 
the logarithmically discretized conduction band. In the 
past few years, the NRG has been extended by divid- 
ing the chain into a system and an environment as in 
the DMRG^; the resulting method, the density-matrix 
numerical renormalization group (DM-NRG), has under- 
gone significant development recentlyj 40-44 While these 
extensions have led to significant improvements, they 
nevertheless have focused primarily on optimizing algo- 
rithms based on a one-dimensional topology. In general, 
the conduction bands are entangled through the impu- 
rity subsystem only. Therefore, it is a natural choice to 
describe the problem as a tree-like structure in which 
the impurity subsystem at the center is surrounded by 
shells of the conduction bath. Such a tree-like topology 
has been utilized to treat the quantum impurity problem 
that occur with dynamic mean field theory (DMFT) i 45 i 46 

The first attempts to apply matrix-product states on 
tree networks were carried out by Otsuka 4 ^ and by Fried- 



man et. al4s, who calculated ground-state properties 
of the spin-1/2 XXZ and Heisenberg chains using the 
DMRG. Subsequently, Lepetit et. al4^ studied the half- 
filled Hubbard model on a Bethe lattice (also known as a 
Cayley tree). All of this work uses a Bethe lattice, whose 
characteristics are that the number of nearest neighbors 
at each node is z, i.e., the coordination number, and that 
closed loops do not occur. Since there is only one path 
between any pair of sites, a DMRG-bascd solution of the 
problem is possible. In this approach, however, z systems 
blocks must be renormalized to a single block at each iter- 
ation step, leading to computational cost which increases 
exponentially with z, hindering systematic DMRG stud- 
ies for large systems and/or large z\ up to now only z = 3 
has been treated. In contrast, the tree tensor network 
(TTNS) methods we introduce here have a much lower 
computational cost because the topology consists of a 
single site and z blocks, whereas the superblock config- 
uration included 2z blocks or z + 1 in previous DMRG 
attempts. In this work, we present calculations on small 
systems in order to demonstrate the viability of the tree 
tensor network method and to compare its accuracy and 
efficiency with existing DMRG methods. 

Another approach to take advantage of the benefits 
of the tree network was formulated in Ref. [5(J . In this 
work, the tree tensor network is formed by placing phys- 
ical sites on the boundary sites only. The remaining in- 
terior sites are virtual; they are only used to transfer 
entanglement up the tree. This tree tensor product state 
is designed to treat models in which sites contribute the 
same amount of entanglement, i.e, have the same value 
of single-site entropy. 

In this work, we form a tree tensor network in which 
all sites in the tree represent physical sites and in which 
entanglement is transferred via the virtual bonds that 
connect the sites. Our motivation is to treat models in 
which physical sites have varying degrees of entangle- 
ment; positions closer to the center of the tree should 
be better suited to represent more entangled sites. An 
additional motivation is to take advantage of the prop- 
erty of the tree tensor network ansatz that the long-range 
correlations differ from the mean-field value polynomially 
with distance rather than exponentially with distance as 
for MPS. In our algorithmic approach to optimize the 
tree tensor network, we use tools similar to those used in 
Refs. [HI and [5(J , but optimize the network site- by-site 
as in the DMRG instead of performing an imaginary time 
evolution. In addition, we explicitly describe how to deal 
with fermions and long-range interactions. 

The paper is organized as follows. In Sec. [Til we de- 
scribe the theoretical background for the TTNS algo- 
rithm and the orbital optimization used during the it- 
erative procedure. Sec. |Hl] is devoted to the analysis 
of the error in the ground-state energy for various spin 
and fermion models as a function of bond dimension for 
TTNS with different coordination numbers. These in- 
clude the 2-d Heisenberg model, the 2-d spinless fermion 
model, the 1-d spinless fermion model in momentum 
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space, and the Beryllium atom as an application in quan- 
tum chemistry. In the latter case, we compare the results 
to those from DMRG calculations. We conclude and dis- 
cuss future prospects in Sec. [IVJ 



II. THEORETICAL BACKGROUND OF THE 
TREE-NETWORK APPROACH 

We approach the problem of finding the ground state of 
strongly correlated systems with long-range interactions 
using a TTNS ansatz of the form 

*)= Ck 1 ...k M \k 1 ,...,k M }- 

ki,...,kM 

Here the coefficients C^...^ M describe a tree tensor net- 
work, i.e., they emerge from contractions of a set of ten- 
sors {Ai, . . . , Am} according to a tree network, as shown 
in Fig. [TJ We associate a tensor with z + 1 indices, 

i^-m]a 1 ...a z ' 

with each vertex m of the network, that is, each tensor 
has z virtual indices a± . . . a z of dimension D and one 
physical index k of dimension d, with z being the coordi- 
nation number of that site. The coefficients Cfc 1 ... / t M are 
obtained by contracting the virtual indices of the tensors 
according to the scheme of a tree tensor network (see 
Fig. [TJ. The structure of the network can be arbitrary 
and the coordination number can vary from site to site. 
The only condition is that the network is bipartite, i.e., 
by cutting one bond, the network separates into two dis- 
joint parts. In the special case z = 2, the one-dimensional 
MPS-ansatz used in DMRG is recovered. 

In a tensor network, entanglement is transferred via 
the virtual bonds that connect the sites. Thus, it is 
preferable to put strongly correlated sites close together, 
i.e., to minimize the number of bonds between them. Evi- 
dently, the diversity of arranging the sites in the network 
increases drastically with increasing coordination num- 
ber z. Also, with a coordination number z > 2 the num- 
ber of virtual bonds required to connect two arbitrary 
sites scales logarithmically with the number of sites M, 
whereas the scaling is linear in M for z = 2 £2- This can be 
seen by considering a Cay ley-tree of depth A, as shown 
in Fig. [TJ The number of sites in the tree is 

M = i+,x>-ir* = ^-^- 2 

and thus, the maximal distance between two sites, 2A, 
scales logarithmically with M for z > 2. This logarith- 
mic scaling is fundamental because, with a MPS ansatz, 
the expectation value of a long-range correlation differs 
from the mean-field result only by a quantity that decays 
exponentially with the distance: 

(t„t„ + a) - (t„)(t„ +a ) oc c" a 



With a TTNS ansatz, the logarithmic scaling counter- 
acts this exponential decay, so that the difference from 
the mean-field result only scales polynomially with the 
distance. 

The way to arrange the physical sites on the network 
is determined by the choice of the basis | fei, &m )• Ob- 
viously, the precision of the ansatz depends critically on 
the choice of the basis. For example, the noninteracting 
Fermi gas has a ground state that is a direct product in 
the momentum space representation. Such a state cor- 
responds to a tree tensor network state with virtual di- 
mension D = 1 at all bonds. In position representation, 
however, a bond dimension that increases exponentially 
with the number of sites would be required. Thus, it will 
be favorable to optimize not only over the tensors in the 
coefficients C7fe 1 ...fc M , but also over the basis \k\, A;m )■ 

In second quantization, | k\, /cm ) denotes the basis 
in occupation number representation, 

\k 1 ,...,ku} = (a\) kl — (J M ) kM |0) 

and a basis transformation is obtained from the canonical 
transformation 

M 

b](U) = ]T U jr at 

r=l 

defined by the M xM unitary matrix U. The whole vari- 
ational ansatz then depends on two sets of parameters: 
the tensors {Ai, . . . , Am} and the unitary U. The goal of 
the algorithm is to find the minimium of the energy with 
respect to these two sets of parameters, i.e., to calculate 

E= min (¥(A 1 ,...,A M ,U)\H\*(A 1 ,...,A M ,U)). 

Ai...A M U 

The idea is to perform the optimizations over the param- 
eter sets {Ai, . . . , Am} and U consecutively and repeat- 
edly until convergence is reached. 

In the following two sections, we describe the two op- 
timization procedures in more detail. 

A. Network Optimization 

First, let us sketch how to optimize the tensors 
{Ai, . . . , Am} while keeping the basis fixed. The mini- 
mization of the energy with the constraint that the norm 
remains constant is equivalent to minimizing the func- 
tional 

F = (* \H\y) - | *) - 1). 

This functional is non-convex with respect to all param- 
eters {Ai, . . . , Am}- However, due to the tensor network 
structure of the ansatz, it is quadratic in the parameters 
A m associated with one lattice site m. Because of this, 
the optimal parameters A m can simply be found by solv- 
ing a generalized eigenvalue problem T-L m A m = EJ\f m A m . 




FIG. 2. (Color online) (a) Separation of the state into z 
blocks plus the site under optimization, as described by 
Eq. <[3j . (b) Natural freedom in the tensor network: insertion 
of a matrix w and w' 1 at one bond leaves the state invariant. 
The contraction of A with w forms the new tensor A' on the 
left hand side; the contraction of B with w -1 forms the new 
tensor B' at the right hand side. 

For a bipartite network, it is always possible to assume 
a gauge condition so that Af m = I, and thus reduce the 
generalized eigenvalue problem to an ordinary one^. We 
will discuss this in more detail later in this section. The 
concept of the algorithm is to do this one-site optimiza- 
tion site-by-site until convergence is reached. 

The challenge that remains is to calculate the effective 
Hamiltonian % m of the eigenvalue problem. In principle, 
this is done by contracting all indices in the expression 
for the expectation value ( \& \H\ \& } except those that 
connect to A m . By interpreting the tensor A m as a dD z - 
dimensional vector A m , this expression can be written 
as 

{*\H\H>)=Al n H m A m . (1) 

Since 

(*\*) = AlAf m A m (2) 
and M m = I, the functional F attains its minimum when 

^m^m Fj Am . 

Due to the bipartite structure of the tensor network, the 
calculation of 7i m can be performed efficiently, i.e., on a 
time that scales polynomially with M and D. Assuming 
that the coordination number z is the same everywhere, 
the scaling will be MdD z+1 . 

This procedure is similar to a DMRG calculation 
with z blocks instead of two, where a block consists of 
all of the sites within one of the branches emerging from 
site m (see Fig. H^a)). The wave function is then formed 



(a)< V |h r iv> = 




FIG. 3. (Color online) (a) Formation of the effective Hamilto- 
nian rj m , r = f)m r <S> hm, r ® f)m,r with respect to the interaction 
h r = r 1 - 7 ' ® t' 15 '. The sites on which the interaction has 
support are marked in red. Each open (filled) circle in the 
tensor network corresponds to the contraction of the layered 
structure of tensors shown in (b). 

as 

D 

i*>= E i¥w..*.>®iO®---®K.>> ( 3 ) 

ai , . . .,a z — 1 

where | <j> l a ) (a = 1, . . . , D) is the basis in environment 
block i (i = 1, . . . , z) and | <p ai ...a z ) is the state of site m. 
Since M m is obtained from the norm ( ^ \ ^ ) by contract- 
ing all tensors except A m [see Eq. ©], it factorizes into 
a tensor product of z matrices Af ml 

where each matrix M m is formed by taking the overlap 
of the basis states in environment block i: 

Obviously, in order to guarantee that VV m = I, the basis 
in each environment block must be orthonormal. This 
is a similar requirement as in the DMRG. In the tree 
tensor network, this may be achieved by assuming an 
appropriate gauge condition. The purpose of this gauge 
condition is to fix the natural freedom in the tensor net- 
work that a matrix and its inverse can be inserted at any 
bond, with the matrix being contracted with the first 
attached tensor and the inverse being contracted with 
the second attached tensor, leaving the network invari- 
ant (see Fig. [2Kb)). The gauge condition for all sites n 
(n ^ m) that assures that Af m = I is 

X/ [An]a'p 2 ...0; lAJaft,...^ = <W • (4) 
kfa...$ % 

Here we take the index a to be outgoing, i.e., the branch 
attached to this index contains site to, and the indices 
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< \|/ | h r | y > = 




FIG. 4. (Color online) Branches with no support (marked 
by dotted lines and circles), which yield the identity when 
contracted and therefore do not have to be taken into account 
in calculating the effective Hamiltonian. 

/?2 ■ ■ • Pz to be incoming, i.e., the attached branches do 
not contain site m (see Fig. [2ja)). 

A stable way to apply this gauge condition is to "or- 
thonormalize" the tensors from outside to inside. That 
is, starting with the tensors on the leaves of the tree, 
we take into account condition (jU by QR-dccomposing 
[A n ] a , i.e., by writing it as 

[A n ] a = IQnla' Ra'a- 
ot' 

The unitary matrix [Qn]^' is the new "orthonormalized" 
tensor at site n. In order to keep the tensor network 
invariant, R a ' a must be contracted with the tensor on the 
first inner shell that is connected to the tensor at site n. 
This procedure continues iteratively with the tensors on 
the first inner shell, the second inner shell, and so on 
until site m is reached. 

Thus, by assuring that the gauge condition is always 
satisfied in the course of the algorithm, the only term 
that must be calculated is the effective Hamiltonian T-L m . 
This term is obtained, as can be gathered from Eq. ([1}, 
by contracting all tensors except A m in the expectation 
value \H\^>). Since the Hamiltonian is a sum of in- 
teraction terms 

R 

H = xx> 

r=l 

with h r being a tensor product of matrices, e.g., h r = 
t^ 7 ) Cg> t' 15 ) for a two-body interaction acting on sites 7 
and 15, the effective Hamiltonian splits into a sum of 
effective Hamiltonians i) m ,r with respect to a single in- 
teraction term h r : 



<v\K\v> = 




FIG. 5. (Color online) Formation of the effective Hamilto- 
nian \) m ,r = l)m,r ® fym,,r ® t)m,r with respect to the fermionic 
interaction h r = a\ais- The sites on which the interaction 
has support are marked in red. 

Due to the structure ^ of the TTNS, each effective 
Hamiltonian f) m r factorizes into a tensor product of z 
matrices 

f)m,r = f) m , r ® • • • ® f)m,r ! 

where each matrix tf m r corresponds to the matrix el- 
ements of h r with respect to the basis in environment 
block i: 

Graphically, the evaluation of ( <p l a \ h r \ </>^ ) corresponds 
to the contraction of a three-layered tensor network ac- 
cording to the structure of the branch in block z, as de- 
picted in Fig. [3] This network can be contracted effi- 
ciently by starting from the leaves and working in the 
inward direction. The numerical effort for contracting 
one additional site into the network scales as dD z+1 , so 
that the total effort scales as dD z+1 times the number of 
sites in the block. 

Of course, if h r has no support in environment block i, 
t) l m r is equal to the identity because the basis is chosen 
to be orthogonal in each environment block. Because of 
this, the calculation simplifies significantly. For exam- 
ple, each two-site interaction has support in at most two 
blocks. This means that at most two effective Hamil- 
tonians tf m r have to be calculated (for each interaction 
term); all others are equal to the identity. Since the or- 
thonormalization of the state is applied iteratively from 
the leaves inwards to the optimized site m, the calcula- 
tion of t) l m r simplifies substantially, as well. Within each 
block, the contraction of all subbranches on which \) l m r 
has no support automatically yields the identity. Thus, 
for the example of a two-site interaction, it is sufficient 
to take into account the sites on the path connecting 



the two sites (see Fig. [4} . The treatment of long-range 
interactions on a tree is therefore not significantly more 
complicated than the treatment of long-range interaction 
on a chain in the DMRG. 

At first glance, it might seem that the advantage of 
"cropping" all subbranches with no support is lost when 
switching to fermions. This is because fermionic inter- 
actions with local support are turned into interactions 
with nonlocal support in the spin picture via the Jordan- 
Wigner transformation. The two-particle fermionic in- 
teraction a^ais, for example, is turned into the spin in- 
teraction 



a 7 ai5 



r < 7 > 



n 

,7<n<15 



z (n) 



r (15) 



that has support on 9 sites, where Z = —a z and rr + , cr_, 
a z denote the Pauli operators (see Fig. [SJ. However, as 
we will show, local fermionic interactions can be treated 
in the tree network with the same effort as local spin in- 
teractions by including the ^-symmetry in the ansatz. 
This is the fermion number parity in the fermionic lan- 
guage. 

The /^-symmetry can be incorporated into the ansatz 
by making the tensors [A^j^ a block-diagonal. That 
is, each virtual index is split into an index pair (ai,Pi), 
with pi carrying the parity information. By stipulating 
that pi © • • • © p z © k = 1 (where © denotes summation 
modulo 2), it is possible to "move" the operator Z that 
acts on the physical index of A m to the virtual parity 
indices pf. 



^kk [Am] aiPl . 



[A Tl 



Z n 



From this relation, it immediatly follows that 
Z®...®Z|#) = |tf) 



(5) 



(6) 



and thus that the state is ^-symmetric. This is because 
the left-hand side of Eq. © corresponds to the contrac- 
tion of Z's to the physical indices of all tensors A m . After 
moving all Z's to the virtual bonds, the operator Z ap- 
pears twice on each bond . Thus, since Z 2 = I, the state 
I \& ) remains unchanged. 

Using the same idea, we can immediately enforce a 
more restrictive symmetry that is fulfilled by all fermionic 
Hamiltonians: the charge symmetry U(l), i.e. the con- 
servation of the number of particles. For this, the tree 
graph has to be made directed (see Fig. [5]), so that all 
sites (except one) have z — 1 incoming and one outgo- 
ing bond. The exception is the sink site with z incom- 
ing bonds. Thus, each virtual index of a tensor A m is 
equipped with the additional information of whether it 
is "incoming" or "outgoing" . We assume that the index 1 
is always the outgoing index in the following. As before, 
each virtual index is split into an index pair (o<j,ni). In 
addition, we require that n\ = 712 + . . . + n z + k. Thus, 
for i = 2, . . . , z, Hi counts the number of particles within 



(a) 




0*1, Hi ) 



(b) (AT 

(a 3 ,n 3 )^-^(a 2 ,n 2 ) 



FIG. 6. (Color online) (a) Possible choice of an ordered tree 
graph with each site having two incoming and one outgoing 
bond. The sink site is marked in black. The tensors A m 
associated with each site with virtual indices consisting of 
pairs (ai,rii) have the structure shown in (b). 



the branch attached to index i. The index m, on the 
other hand, is equal to the number of particles in all the 
branches plus the number of particles at site m. Since 
the parity can be immediately derived from the particle- 
number information, a similar relation as Eq. ([3]) holds, 
namely 



^kk [^Lmr 



[A 



k 



(7) 

with Znn = <5fm(— !)"■ Thus, as before, Z acting on 
the physical index can be moved to the virtual particle- 
number bonds and, since Z 2 = I, the state is Z2- 
symmetric. 

The advantage of taking into account the particle num- 
ber is twofold: on the one hand, the block structure of 
the tensors reduces the numerical effort considerably: a 
bond-dimension of \ — ND can be treated with an effort 
of order N Z D Z+1 instead of x 2+1 for a non-symmetric 
ansatz. Here % is the full bond dimension, so it is equiv- 
alent to the number of block states kept in a DMRG 
procedure. On the other hand, the calculation of the 
effective Hamiltonians f) mjJ , with respect to an inter- 
action with support on a few sites only, e.g., a two- 
particle or four-particle interaction, is simplified. The 
main idea is depicted in Fig. [7] for the interaction a\a\^: 
with an appropriately chosen numbering of the fermions, 
each subbranch that has no fermionic support cither has 
only identities acting on the sites or only operators Z . 
The subbranches including only identities simplify to the 
identity because of the orthonormalization of the state. 
The Z operators can be moved, using relation (|7|), to 
the virtual bonds, and all of them except one cancel (see 
Figs. O and [7]). What remains is a subbranch that in- 
cludes only identities, which reduces to the identity be- 
cause of the orthonormalization of the state. Thus, for a 
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< \|/ | h r | y > = 




\ / N \ ' \12 
\20 0--~019 \ / 013/ 

\ 621 \, 614 / 



FIG. 7. (Color online) Formation of the effective Hamiltonian 
flm,r = f)m, r ® f)m,r ® f)m, r with respect to the fermionic inter- 
action ft, r = a\aiB with parity symmetry taken into account. 
The sites on which the interaction has support are marked in 
red. All branches marked by dotted lines and circles yield the 
identity when contracted. The parity operator Z is contracted 
to the virtual bond. 



fermionic two-site interaction, it is sufficient to take into 
account the path connecting the two sites. In this way, 
the treatment of long-range fermionic interactions is fea- 
sible with the same numerical effort as the treatment of 
long-range spin interactions. 

Using the same scheme as for minimizing the energy, 
the efficient simulation of time evolution is also possible. 
For this, the time evolution is split up into small steps 
of duration 5t. For 5t ^ 1 and a starting state | <3/o } of 
the form of a TTNS, the time-evolved state e~ iHSt \ * ) 
is a TTNS as well. However, the virtual dimension is 
multiplied by a factor k > 1. In order to prevent the 
bond dimension from increasing exponentially with time, 
the TTNS has to be approximated by a TTNS with a 
reduced virtual dimension | ) after every time step, i.e., 
the functional 

K= \\e-' lHt \^ ) - |*)|| 2 

must be minimized. The optimization with respect to 
a single site m leads to the system of linear equations 
A/" m j4 m = w m (where Af m can again be set equal to iden- 
tity by using the appropriate gauge) with 

(*|e- lHt |* > =Alw m 

and thus can be performed efficiently with the scaling 
MD Z+1 , as before. Clearly, all previous considerations 
regarding particle-number symmetry can also be adapted 
to the case of time evolution. 



B. Orbital Optimization 

As mentioned previously, the entanglement properties 
of the system depend critically on the choice of the basis. 
Our goal is to find a basis in which entanglement is local- 
ized as much as possible at the sites of the tree network. 
Such a choice of basis guarantees that a given precision 
can be attained with a smaller virtual dimension D, and 
thus with less computational effort. In QC-DMRG the 
optimization of the basis is fundamental and has been 
used in much previous worki 32 i 51 ~ — 

In contrast to previous work, our approach aims to find 
the optimal basis that can be obtained by a canonical 
transformation of the fermionic modes. The canonical 
transformation is defined by a M x M unitary matrix 
U . Since the number of parameters is relatively small, 
a gradient search applied to the expectation value of the 
energy, 

E(U) = (*(U)\H\*(U)), 

is certainly feasible. Since E{U) is a non-convex function 
of the parameters U, it is a highly non-trivial problem to 
find the absolute minimum. The idea is to perform the 
gradient search multiple times in the course of the algo- 
rithm, e.g., after each optimization sweep of the tensor 
network, and improve the energy at each gradient search 
by only a small amount. This will eventually adapt the 
orbitals optimally to the tree tensor network. 

There are two ways to implement the basis transfor- 
mation: one based on the state and the other based on 
the Hamiltonian. We have applied the basis transforma- 
tion to the Hamiltonian. For the fermionic Hamiltonian 
with long-range interactions, 

H = J2 T v a l a 3 + Vijkl4 af j a kai , 

ij ijkl 

that appears, e.g., in quantum chemistry, in momentum- 
space descriptions of the Hubbard model, or in descrip- 
tions of the Hubbard model on higher dimensional lat- 
tices, the function E(U) can be be expressed as 

E(U) = ]T f fflMaj) + £ V{U) m {a\a]a kai ) 

ij ijkl 

with 

f(U) = UTU f 

V(U) = {U ®U)V{U ®U)\ 

The correlation functions (aj a j) and (ala^a^ai) are cal- 
culated with respect to the original state and are not 
dependent on the parameters in U. With the function 
E{U) in this form, its gradient can be calculated explic- 
itly. Both quantities can be evaluated efficiently for dif- 
ferent parameter sets U, which makes the gradient search 
feasible. 

The gradient search that is used for the orbital opti- 
mization typically finds the local minima in the vicinity 
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FIG. 8. (Color online) Relative error in the energy for Heisen- 
berg model on a 4 x 4-lattice (a) and a 6 x 6-lattice (b) as a 
function of the optimization steps. The calculations were per- 
formed with a fixed virtual dimension of D — 4 and tree ten- 
sor networks with coordination numbers 2 = 2 (blue curve), 
z — 3 (green curve) and 2 = 4 (red curve). 

of the starting point of the optimization. In order to 
avoid falling into the same minimum, we shift the start- 
ing point in a random direction by an appropriately cho- 
sen amount. This has the consequence that the gradient 
search generally falls into another local minimum that 
may have a higher energy. However, it assures that the 
orbitals change considerably and that, since the orbital 
optimization is performed repeatedly interspersed with 
network optimization, the energy decreases notably in 
the course of the algorithm. 

III. NUMERICAL RESULTS 

We have applied the algorithm consisting of two op- 
timization procedures described above, the optimization 
of the tensor network and the optimization of the ba- 
sis according to this network, to several models. In this 
section, we discuss the results. 

A. 2-d Heisenberg model 

First, let us consider the tree tensor network optimiza- 
tion only and show using the 2-d Heisenberg model how 
results improve with a TTNS ansatz compared to a one- 
dimensional MPS ansatz. The relative error in the energy 
of a system consisting of 4 x 4 spins as a function of the 
optimization step is shown in the left hand side of Fig. [8] 
for the MPS ansatz and the TTNS ansatz with z = 3 
and z = 4. In all calculations, a fixed virtual dimension 
of D = 4 is used. We have assigned a spin to each node in 



FIG. 9. (Color online) Relative error in the energy for the 
interacting spinless fermion model on (a) a 4 x 4-lattice and 
(b) a 6 x 6-lattice as a function of the optimization steps. The 
calculations were performed with a fixed virtual dimension of 
D = 4 and tree tensor networks with coordination numbers 
2 = 2 (blue curve) and 2 = 3 (green curve). The chosen 
parameters are J — 1, U = 0.5 and the number of fermions is 
fixed to N = 3. 

the tree network in such a way that two arbitrary spins 
are connected by the smallest possible number of bonds. 
As can be seen, the precision increases considerably with 
increasing coordination number. 

For larger systems such as for 6x6 spins, shown in 
the right hand side of Fig. [51 the energy is plotted as 
a function of the optimization steps for z = 2, z = 3, 
and 2 = 4. The virtual dimension is fixed to D = 4. As 
before, the energy decreases considerably and approaches 
the PEPS result for D = 4, indicated by the dotted line 
which, however, is obtained with a much larger numerical 
efforts 



B. 2-d interacting spinless fermions 

The TTNS Ansatz also perfoms well on two- 
dimensional fermionic models. We have applied the 
TTNS algorithm to a system of interacting fermions on 
a two-dimensional lattice, described by the Hamiltonian 

H = - J 4 a J + U X] ™ lfl i 
<i,j> <i,j> 

with hi = a\a,i. The boundary conditions are assumed to 
be periodic. As can be gathered from Fig. |H1 the ground- 
state energy improves with increasing coordination num- 
ber z. The figure shows the relative error in the energy as 
a function of the iteration step for z = 2 and z = 3 for a 
4 x 4-latticc in the left part, and the ground state energy 
as a function of the iteration step for a 6 x 6-lattice in the 
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FIG. 10. (Color online) Relative error in the energy for in- 
teracting spinless fermions with M = 7 sites as a function of 
the optimization steps for the TTNS ansatz with z = 2 (a) 
and 2 = 3 (b). Here U = 1, the particle number is fixed to 
TV = 3, and the virtual dimension is assumed to be D — 2. 
Results with orbital optimization are compared to the results 
obtained in the real-space basis and in the momentum-space 
basis. The steps at which orbital optimizations are performed 
are marked with arrows. 



right part. For the calculations, we have used the virtual 
dimension D = 4. We have chosen the parameters J = 1 
and U = 0.5 and have fixed the number of fermions to 
TV = 3. 

Thus, a TTNS Ansatz might be useful for the study 
of higher dimensional models of small size because the 
effective long-range interactions are represented better in 
a tree than in a chain and the numerical effort is relatively 
low compared to a PEPS calculation. 



C. 1-d interacting spinless fermions 

In order to assess the effectiveness of the orbital op- 
timization, we have applied the algorithm to a simple 
fcrmionic model, the one-dimensional interacting spinless 
fcrmion model 



H 



M 

E 



M 



a] a*. 



h.c. 



where hi = a\a,i. We assume periodic boundary condi- 
tions, i.e., om+i = a>i- This model can be mapped to the 
XXZ spin chain via a Jordan- Wigner transformation. In 
this model, it is known that the choice of the basis has 
a big effect on the precision of the DMRG or TTNS cal- 
culation. For U — > oo, the ground state is a product 
state in the position representation, and thus optimally 
represented by a tensor network with D = 1 in this ba- 
sis. For U — 0, on the other hand, the ground state can 



FIG. 11. (Color online) The relative error in the ground- 
state energy for the Beryllium atom as a function of DMRG 
iteration steps for various values of the number of DMRG 
block states. The dashed line corresponds to the Hartree- 
Fock energy. 



be represented as a direct product in momentum space. 
Thus, the momentum-space basis is clearly best suited in 
this limit. For intermediate U, one might expect a basis 
intermediate between real and momentum space to rep- 
resent the entanglement properties optimally; our goal is 
to find such a basis automatically by carrying out orbital 
optimization. 

The results of calculations incorporating the optimiza- 
tion procedure are displayed in Fig. [TU] for z = 2 and 
z = 3, both on systems of M = 7 sites and N = 3 
particles with a fixed virtual dimension of D = 2 and 
parameters U = 1. As before, the relative error in the 
energy as a function of the optimization step is plotted. 
For comparison, calculations performed in the position 
representation and in the momentum representation are 
also shown. As can be seen, the energy improves signifi- 
cantly in the course of the optimization; the improvement 
is three orders of magnitude for the z = 3 case. 



D. Quantum Chemical systems 

In this section, we compare numerical results for quan- 
tum chemical systems obtained using the QC-DMRG 
and TTNS methods. In these applications, the electron- 
electron correlation is taken into account by an iterative 
procedure that minimizes the Rayleigh quotient corre- 
sponding to the Hamiltonian describing the electronic 
structure of the molecule, given by 



ijklaa' 



aic 



(8) 
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FIG. 12. (Color online) Relative error in the ground-state en- 
ergy for the Be atom as a function of the optimization step for 
D = 2. The results obtained using a TTNS ansatz with z = 3 
plus orbital optimization are compared to the TTNS-results 
without orbital optimization For comparison, the results with 
a MPS ansatz (corresponding to z = 2) are also included. The 
steps at which orbital optimizations are performed are marked 
with arrows. 

and thus determines the full-CI wave function. In 
Eq. ([8]), Tij denotes the matrix elements of the one- 
particle Hamiltonian, which is comprised of the kinetic 
energy and the external electric field of the nuclei, and 
Vijki stands for the matrix elements of the electron re- 
pulsion operator, defined as 

V ijkl = [ d 3 x 1 d 3 x 2 <S>*{x 1 )m(x 2 ) z; 1 z; $ k (x2)$l(xi) ■ 

J J Xi~X 2 

We obtain the Hartree-Fock orbitals in a given basis of 
Gaussian orbitals and transform the matrix elements TJy 
and V^ki to the Hartree-Fock basis using the MOLPRO 
program package,— which we also use to obtain the full- 
CI energies used as a benchmark. 

In the QC-DMRG, a one-dimensional chain is built up 
from the atomic or molecular orbitals obtained from a 
suitable mean-field or MCSCF calculation. The tree net- 
work is constructed similarly, but there is greater free- 
dom to form the proper structure of the network. The 
two-orbital mutual informationiS provides a good start- 
ing configuration. A general approach to reduce entan- 
glement is to form the network by placing the highly 
correlated orbitals at or near the center of the tree and 
less correlated orbitals at the boundary sites. 

In Fig.QTJ we plot the relative error in the ground-state 
energy for the Beryllium atom as a function of DMRG it- 
eration steps for various fixed values of the DMRG block 
states. Corresponding data gathered after the fourth 
DMRG sweep is shown in Fig. Q2J In this calculation, 
four electrons have been placed on eight orbitals. The 



FIG. 13. (Color online) Relative error in the ground-state 
energy for the Beryllium atom as a function of the full bond 
dimension \ obtained using (a) the DMRG and (b) as a func- 
tion of the virtual dimension D obtained from a TTNS with 
z — 2 and 2 = 3. 

Hartree-Fock energy is -14.351880250000, while the full- 
CI energy is -14.37016558404629. Figure. ED depicts the 
relative error as a function of the optimization step for 
the Be atom, calculated using a TTNS ansatz with co- 
ordination number z = 3. The cases with and without 
orbital optimization are considered separately, and the 
z = 2 results are included for comparison. As can be 
seen, the relative error is considerably smaller for z = 3 
than for z = 2, and there is a significant gain in preci- 
sion when orbital optimization is taken into account. In 
Fig. [T31 we display the dependence of the z = 2 and 3 
calculations on bond dimension D. As can be seen, the 
DMRG and TTNS calculations with z = 2 yield similar 
accuracy, while the z — 3 TTNS calculation is signifi- 
cantly more accurate for a given bond dimension. 

It should be noted that there is currently a discrepancy 
between the QC-DMRG and the TTNS calculations in 
the speed of convergence with optimization step. Since 
our QC-DMRG code is highly optimized, the calculation 
converges faster than that in the present version of the 
TTNS method. As can be seen in Fig. [TTJ for example, 
the QC-DMRG ground-state energies are always below 
the Hartree-Fock energy, and the active space is extended 
dynamically using Cl-based expansion techniques (CI- 
DEAS) 5 (DEAS). In general, there is no fundamental 
difficulty in incorporating optimizations such as the CI- 
DEAS into the TTNS method. 



IV. CONCLUSION 

In this paper, we have described and applied a method, 
the tree tensor network state method, to treat strongly 
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correlated systems with long-range interactions on a tree 
network with arbitrary coordination number z. Our ap- 
proach is based on a tensor product state ansatz that 
generalizes a DMRG-like matrix product state to z rather 
than two blocks. The number of virtual bonds required 
to connect two arbitrary sites scales logarithmically with 
the number of sites in TTNS, in contrast to the linear 
scaling of a one-dimensional topology^ In this sense, our 
TTNS method has a lower computational cost than cur- 
rently used DMRG-based methods. We have also incor- 
porated optimization of the single-particle orbitals in our 
method, treating the case of a procedure for tranforming 
the basis that smoothly interpolates between real space 
and momentum space. We have tested our method us- 
ing numerical calculations on various systems with local 
and nonlocal interactions, including the two-dimensional 
Heisenbcrg lattice, the momentum-space version of the 1- 
d interacting spinless fermion model, and small quantum 
chemical systems. For the quantem chemical systems, we 
have compared TTNS results to those of DMRG calcu- 
lations. 

Since the TTNS approach is defined on a bipartite net- 
work, previous algorithmic developments and optimiza- 
tions procedures developed in the context of the quan- 
tum chemistry version of DMRG can also be integrated 



into the TTNS method. Such optimizations include dy- 
namic adjustment of the bond dimension^ 5 - orbital opti- 
mization, 32 i 51 - 54 i 56 i 57 an d initialization procedures based 
on CI expansions ; 11 ! 14 ' 58 among others. Incorporation of 
these elements into the TTNS approach will be carried 
out in future work. In light of the promising features of 
the new method, we expect it to provide a viable alter- 
nate means of treating atoms and molecules efficiently in 
the near future. 
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